# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Finite differences stencil inteface used in various numerical methods.
* :class:`~hysop.numerics.stencil.stencil.Stencil`
* :class:`~hysop.numerics.stencil.stencil_generator.StencilGenerator`
"""
import copy
import fractions
import gzip
import itertools as it
import math
import os
import numpy as np
import scipy as sp
import sympy as sm
try:
import cPickle as pickle
except:
import pickle
from hysop.numerics.stencil.stencil import CenteredStencil, Stencil
from hysop.tools.cache import load_data_from_cache, update_cache
from hysop.tools.io_utils import IO
from hysop.tools.misc import prod
from hysop.tools.numerics import F2Q, MPFR, MPQ, MPZ, mpq, mpqize, mpz
from hysop.tools.sympy_utils import (
build_eqs_from_dicts,
factor_split,
tensor_symbol,
tensor_xreplace,
)
from hysop.tools.htypes import extend_array
try:
import flint
has_flint = True
except ImportError:
flint = None
has_flint = False
[docs]
class StencilGeneratorConfiguration:
_dx = ("dx", "dy", "dz")
def __init__(self):
self.dim = 1
self.dtype = MPQ
self.dx = sm.Symbol("dx")
self.user_eqs = {}
self.derivative = 1
self.order = 1
self.mask_type = StencilGenerator.CROSS
self._mask = None
self._format_inputs()
[docs]
def copy(self):
obj = StencilGeneratorConfiguration()
for var in [
"dim",
"dtype",
"dx",
"user_eqs",
"derivative",
"order",
"mask_type",
"_mask",
]:
setattr(obj, var, copy.deepcopy(getattr(self, var)))
return obj
[docs]
def shape(self):
if (self.derivative is None) or (self.order is None):
raise RuntimeError("First set derivative and order with configure!")
return self.derivative + self.order
[docs]
def L(self, origin):
StencilGeneratorConfiguration._check_origin(origin, self.dim)
shape = self.shape()
return origin
[docs]
def R(self, origin):
StencilGeneratorConfiguration._check_origin(origin, self.dim)
shape = self.shape()
return shape - origin - 1
[docs]
def mask(self, origin, custom_mask=None):
StencilGeneratorConfiguration._check_origin(origin, self.dim)
self._check_and_raise()
if self.mask == StencilGenerator.CUSTOM:
if (custom_mask is None) and (self._mask is None):
raise ValueError("Custom mask should not be None!")
elif custom_mask is not None:
mask = custom_mask.copy()
self._mask = mask
else:
mask = self._mask
else:
mask = StencilGeneratorConfiguration._genmask(
origin, self.shape(), self.mask_type
)
return mask
[docs]
def symbolic_derivatives(self, extra_size=0):
df, df_vars = tensor_symbol(prefix="f", shape=self.shape() + extra_size)
return df, df_vars
[docs]
def symbolic_stencil(self, origin):
StencilGeneratorConfiguration._check_origin(origin, self.dim)
S, svars = tensor_symbol(
prefix="S", mask=self.mask(origin), shape=self.shape(), origin=origin
)
return S, svars
def _format_inputs(self):
if (
isinstance(self.dx, sm.Symbol)
or isinstance(self.dx, sm.Expr)
or np.isscalar(self.dx)
):
self.dx = np.asarray([self.dx] * self.dim)
dim = self.dim
self.derivative = extend_array(self.derivative, dim, dtype=int)
self.order = extend_array(self.order, dim, dtype=int)
self.dx = extend_array(self.dx, dim, dtype=object)
def _check_and_raise(self):
self._format_inputs()
if self.dim < 1:
raise ValueError("Dim < 1!")
if self.dtype not in [MPQ, float, np.float32, np.float64]:
raise TypeError(f"Invalid dtype {self.dtype}!")
if not isinstance(self.dx, np.ndarray):
raise TypeError("Invalid type for dx!")
if self.dx.size != self.dim:
raise ValueError("Invalid size for dx!")
if not isinstance(self.user_eqs, dict):
raise TypeError("Invalid type for user_eqs!")
for k in self.user_eqs:
if not isinstance(k, sm.Symbol):
raise TypeError(f"Invalid type for key {k} in user_eqs!")
if (self.derivative < 0).any():
raise ValueError("derivative < 0!")
if (self.order < 1).any():
raise ValueError("order < 1!")
@staticmethod
def _check_origin(origin, dim):
if not isinstance(origin, np.ndarray):
raise TypeError("Origin is not a np.ndarray!")
if origin.size != dim:
raise ValueError("Bad dimensions for origin!")
if (origin < 0).any():
raise ValueError("Origin component < 0!")
@staticmethod
def _genmask(origin, shape, mask):
dim = origin.size
if mask == StencilGenerator.CROSS:
mask = np.zeros(shape, dtype=bool)
for d in range(dim):
access = [slice(origin[dd], origin[dd] + 1) for dd in range(dim)]
access[d] = slice(None)
access = tuple(access)
mask[access] = True
mask[tuple(access)] = True
elif mask == StencilGenerator.DIAG:
mask = np.ones(shape, dtype=bool)
for d in range(dim):
access = [slice(origin[dd], origin[dd] + 1) for dd in range(dim)]
access[d] = slice(None)
access = tuple(access)
mask[access] = False
mask[tuple(access)] = False
mask[origin] = True
elif mask == StencilGenerator.DENSE:
mask = np.ones(shape, dtype=bool)
else:
raise NotImplementedError("Mask not implemented yet!")
return mask
def __str__(self):
ss = """
StencilGeneratorConfiguration
dim: {}
dtype: {}
dx: {}
user_eqs: {}
derivative: {}
order: {}
mask_type: {}
mask: {}
shape: {}
""".format(
self.dim,
self.dtype,
self.dx,
self.user_eqs,
self.derivative,
self.order,
self.mask_type,
self._mask,
self.shape(),
)
return ss
[docs]
class StencilGenerator:
"""
Generate various stencils.
Results are cached and compressed locally.
"""
# Stencil masks
DENSE = 0
CROSS = 1
DIAG = 2
CUSTOM = 99
[docs]
@classmethod
def cache_file(cls):
_cache_dir = IO.cache_path() + "/numerics"
_cache_file = _cache_dir + "/stencil.pklz"
return _cache_file
def __init__(self, config=None, cache_override=False):
"""
StencilGenerator to generate stencils of dimension :attr:`dim` and type :attr:`dtype`.
Results are cached and compressed locally.
Attributes
----------
dim : int
Dimension of the generated stencils.
dtype : {np.float32, np.float64, :class:`MPQ`}
Datatype of the generated stencils.
Parameters
----------
dim : int
Dimension of the generated stencils.
dtype : {np.float32, np.float64, :class:`MPQ`}
Datatype of the generated stencils.
Other Parameters
----------------
cache_override : bool
If set, overwrite already cached stencils.
Raises
------
ValueError
Raised if dim<1 or if dtype in [np.float32, np.float64] and dim!=1.
TypeError
Raised if dtype is not a valid type.
Notes
-----
Depending on data type :attr:`dtype` a different solver can be used.
For fast generation of 1D stencils, use `np.float32` or `np.float64`.
For n-dimentional exact fractional coefficients stencils use the default `mpq` type.
Exact fractional results (mpq) are cached locally.
Examples
--------
Generate 1D approximative forward first derivative stencil of order 1 and 2.
>>> import numpy as np
>>> from hysop.numerics.stencil.stencil_generator import StencilGenerator
>>> generator = StencilGenerator().configure(dim=1,dtype=np.float64)
>>> s0 = generator.generate_approximative_stencil(origin=0, derivative=1, order=1)
>>> s1 = generator.generate_approximative_stencil(origin=0, derivative=1, order=2)
>>> print('{}\\n{}'.format(s0.coeffs,s1.coeffs))
[-1. 1.]
[-1.5 2. -0.5]
Generate 1D exact backward first derivative stencil of order 1 and 2.
>>> from hysop.numerics.stencil.stencil_generator import StencilGenerator
>>> generator = StencilGenerator().configure(dim=1)
>>> s0 = generator.generate_exact_stencil(origin=-1, derivative=1, order=1)
>>> s1 = generator.generate_exact_stencil(origin=-1, derivative=1, order=2)
>>> print('{}\\n{}'.format(s0.coeffs,s1.coeffs))
[-1 1]
[1/2 -2 3/2]
Generate centered 2D Laplacian of order 2:
>>> import sympy as sm
>>> from hysop.numerics.stencil.stencil_generator import StencilGenerator
>>> generator = StencilGenerator().configure(dim=2, derivative=2)
>>> laplacian = generator.generate_exact_stencil(origin=1, order=2, dx=(sm.Symbol('dx'),)*2)
>>> print(laplacian.coeffs)
[[0 1 0]
[1 -4 1]
[0 1 0]]
Similarly, with different custom dx:
>>> laplacian = generator.generate_exact_stencil(origin=1, order=2)
>>> print(laplacian.coeffs)
[[0 dx**(-2) 0]
[dy**(-2) (-2*dx**2 - 2*dy**2)/(dx**2*dy**2) dy**(-2)]
[0 dx**(-2) 0]]
See Also
--------
:class:`Stencil`
"""
config = StencilGeneratorConfiguration() if (config is None) else config
if not isinstance(config, StencilGeneratorConfiguration):
raise TypeError(
"config is not an instance of StencilGeneratorConfiguration!"
)
self._config = config
self._cache_override = cache_override
@staticmethod
def _hash_algorithm():
import hashlib
return hashlib.new("sha512")
@classmethod
def _hash_equations(cls, eqs):
hasher = cls._hash_algorithm()
for e in eqs:
hasher.update(str(e).encode("utf-8"))
return hasher.hexdigest()
[docs]
def get_config(self):
"""
Get current stencil generator configuration state.
"""
return self._config
[docs]
def cache_override(self):
"""
Get caching override flag.
"""
return self._cache_override
[docs]
def generate_approximative_stencil(self, origin, **kargs):
"""
Generate a stencil using hardware floating point arithmetic (fast).
dtype can be np.float16, np.float32, np.float64, MPQ
"""
config = self._config.copy()
config.configure(**kargs)
dim = config.dim
dtype = config.dtype
if dim != 1:
raise ValueError("Bad dimension for approximation stencil generation!")
if has_flint:
solve_dtype = flint.fmpq
elif dtype not in [np.float16, np.float32, np.float64]:
solve_dtype = np.float64
else:
solve_dtype = dtype
dx = config.dx[0]
k = config.derivative[0]
order = config.order[0]
N = config.shape()[0]
origin = StencilGenerator._format_origin(origin, N)
L = config.L(origin)
R = config.R(origin)
if k == 0:
return Stencil([1], [0], 0, dx=dx, error=None)
A = np.empty((N, N), dtype=solve_dtype)
b = np.empty(N, dtype=solve_dtype)
for i in range(N):
b[i] = solve_dtype(int(i == k))
for j in range(N):
A[i, j] = solve_dtype(int((j - origin) ** i))
try:
if has_flint:
coeffs = A.ravel()
Afmpq = flint.fmpq_mat(*(A.shape + (coeffs,)))
Afmpq_inv = Afmpq.inv()
Ainv = np.asarray(Afmpq_inv.entries()).reshape(A.shape)
S = Ainv.dot(b)
else:
S = sp.linalg.solve(A, b, overwrite_a=True, overwrite_b=True)
except:
print("\nError: Cannot generate stencil (singular system).\n")
raise
S *= math.factorial(k)
actual_dtype = type(S.ravel()[0])
target_dtype = dtype
if actual_dtype != target_dtype:
if target_dtype in [np.float16, np.float32, np.float64]:
if has_flint and (actual_dtype is flint.fmpq):
def convert(x):
return target_dtype(float(int(x.p)) / float(int(x.q)))
S = np.vectorize(convert)(S)
else:
S = S.astype(target_dtype)
elif target_dtype == MPQ:
if has_flint and (actual_dtype is flint.fmpq):
def convert(x):
return mpq(mpz(x.p.str()), mpz(x.q.str()))
else:
def convert(x):
frac = fractions.Fraction(x).limit_denominator((1 << 32) - 1)
return mpq(frac.numerator, frac.denominator)
S = np.vectorize(convert)(S)
else:
RuntimeError("Type conversion not implemented yet.")
return Stencil(S, origin, order, factor=1 / (dx**k), dx=dx)
[docs]
def generate_exact_stencil(self, origin, **kargs):
"""
Generate a stencil using quotients (slower).
Yield exact solutions by using symbolic calculation
and quotient arithmetic.
Do not use on high order stencils.
"""
config = self._config.copy()
config.configure(**kargs)
config._check_and_raise()
dx = config.dx
dim = config.dim
dtype = config.dtype
order = config.order
derivative = config.derivative
N = config.shape()
origin = StencilGenerator._format_origin(origin, N)
L = config.L(origin)
R = config.R(origin)
df, df_vars = config.symbolic_derivatives(4)
S, svars = config.symbolic_stencil(origin)
user_eqs = config.user_eqs
if all(derivative == 0):
return Stencil([1], [0], 0, dx=dx, error=None)
if len(user_eqs) == 0:
for i, d in enumerate(derivative):
if dim > 1:
access = [slice(0, 1) for _ in range(dim)]
access[i] = slice(d, d + 1)
access = tuple(access)
user_eqs[df[access].ravel()[0]] = 1
else:
user_eqs[df[d]] = 1
def taylor(df, dx, N):
expr = 0
for n in range(max(N)):
expr += taylorn(df, dx, n, N)
return expr
def taylorn(df, dx, n, N):
dim = dx.size
def preficate(it):
return (it <= N).all() and sum(it) == n
nn = range(n + 1)
itd = it.product(nn, repeat=dim)
itd = filter(preficate, itd)
expr = 0
for der in itd:
expr += taylorn_term(df, dx, der)
return expr
def taylorn_term(df, dx, der):
fact = np.asarray([math.factorial(d) for d in der])
return df[der] * prod((dx**der) / fact)
expr = 0
for ids in np.ndindex(*N):
offset = ids - origin
expr += S[ids] * taylor(df, offset * dx, N)
fact = factor_split(expr, df_vars)
eqs = build_eqs_from_dicts(fact, user_eqs)
svars = list([_ for _ in svars if not (_ == 0)])
stencil = Stencil(S, origin, order)
i = 0
sol = {}
cache_file = self.cache_file()
while stencil.variables().intersection(svars):
eqs_key = StencilGenerator._hash_equations(eqs)
nsol = load_data_from_cache(cache_file, eqs_key)
if (self._cache_override) or (nsol is None):
nsol = sm.solve(eqs, svars)
update_cache(cache_file, eqs_key, nsol)
if not nsol:
break
sol.update(nsol)
S = tensor_xreplace(S, sol)
err = 0
while err == 0:
for ids in np.ndindex(*N):
offset = ids - origin
err += S[ids] * taylorn(df, offset * dx, max(N) + i, N + i + 1)
if err != 0:
order += 1
i += 1
eqs = []
for k, eq in factor_split(err, df_vars).items():
if isinstance(eq, int):
continue
eq = sm.simplify(eq.xreplace(sol))
if eq.atoms(sm.Symbol).intersection(svars):
eqs.append(eq)
stencil = Stencil(S, origin, order, dx=dx, error=err)
if (dx == dx[0]).all() and dx[0] != 1:
stencil.refactor(dx[0] ** (-derivative[0]))
return stencil
@staticmethod
def _format_origin(origin, shape):
dim = shape.size
origin = (
np.asarray([origin] * dim, dtype=int)
if np.isscalar(origin)
else np.asarray(origin)
)
if origin.size != dim:
raise ValueError("Bad input dimensions!")
return (origin + shape) % shape
[docs]
class CenteredStencilGenerator(StencilGenerator):
"""
Generate various centered stencils.
Results are cached and compressed locally.
"""
def __init__(self, config=None, cache_override=False):
super().__init__(config, cache_override)
[docs]
def generate_exact_stencil(self, **kargs):
config = self._config.copy()
config.configure(**kargs)
shape = config.shape()
origin = (shape - 1) // 2
stencil = super().generate_exact_stencil(origin=origin, **kargs)
if stencil.is_centered():
return CenteredStencil.from_stencil(stencil)
else:
raise RuntimeError("Generated stencil is not centered!")
[docs]
def generate_approximative_stencil(self, **kargs):
config = self._config.copy()
config.configure(**kargs)
shape = config.shape()
origin = (shape - 1) // 2
stencil = super().generate_approximative_stencil(origin, **kargs)
if stencil.is_centered():
return CenteredStencil.from_stencil(stencil)
else:
raise RuntimeError(f"Generated stencil is not centered: {stencil.coeffs}")
if __name__ == "__main__":
from hysop.tools.contexts import printoptions
sg = StencilGenerator()
with printoptions(precision=4):
sg.configure(dim=1, derivative=0, order=2, dtype=np.float64)
print("\ndim=1, 0th derivative, np.float64, approximative, shape=(5,):")
for i in range(sg.get_config().shape()[0]):
stencil = sg.generate_approximative_stencil(origin=i).reshape((5,))
print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}")
sg.configure(dim=1, derivative=1, order=2, dtype=np.float64)
print('\ndim=1, 1st order first derivative, np.float64, approximative, shape=(5,):')
for i in range(sg.get_config().shape()[0]):
stencil = sg.generate_approximative_stencil(origin=i).reshape((5,))
print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}")
sg.configure(dim=1, derivative=2, order=2, dtype=np.float64)
print('\ndim=1, 2nd order first derivative, np.float64, approximative:')
for i in range(sg.get_config().shape()[0]):
stencil = sg.generate_approximative_stencil(origin=i)
print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}")
print("\ndim=1, 2nd order first derivative, exact:")
sg.configure(dtype=MPQ, derivative=1)
for i in range(sg.get_config().shape()[0]):
stencil = sg.generate_exact_stencil(origin=i)
print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}")
print("\ndim=1, 2nd order second derivative, exact:")
sg.configure(dtype=MPQ, derivative=2)
for i in range(sg.get_config().shape()[0]):
stencil = sg.generate_exact_stencil(origin=i)
print(f" origin: {i} => {stencil.factor} . {stencil.coeffs}")
print("\n 2D Laplacian, 2nd order")
sg.configure(dim=2, order=2)
laplacian = sg.generate_exact_stencil(origin=1)
print(laplacian.coeffs)
laplacian = sg.generate_exact_stencil(
origin=1, dx=[sm.Symbol("dy"), sm.Symbol("dx")]
)
print("\n", laplacian.coeffs)
df, _ = sg.get_config().symbolic_derivatives()
user_eqs = {df[0][2]: sm.Symbol("Cx"), df[2][0]: sm.Symbol("Cy")}
stencil = sg.generate_exact_stencil(origin=1, user_eqs=user_eqs)
print("\n", stencil.coeffs)
sg = CenteredStencilGenerator()
sg.configure(derivative=2)
print("\nCentered second order derivative stencils:")
for i in range(1, 4):
stencil = sg.generate_approximative_stencil(order=2*i, dtype=np.float16)
print(f' {stencil.coeffs}')
print()
for i in range(1, 4):
stencil = sg.generate_approximative_stencil(order=2 * i, dtype=MPQ)
print(f" {stencil.coeffs}")
print()
for i in range(1, 4):
stencil = sg.generate_exact_stencil(order=2*i)
print(f' {stencil.coeffs}')